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Abstract. We numerically study a two-dimensional granular gas of rigid disks where 
an external driving force is applied to each particle in such a way that the system 
is driven into a steady state by balancing the energy input and the dissipation due 
to inelastic collisions. Two intruder particles embedded in this correlated medium 
experience a fluctuation-induced force - that is itself a fluctuating quantity - due 
to the confinement of the hydrodynamic fluctuations between them. We find that 
the probability distribution of this force is a Gaussian centered on a value that is 
proportional to the steady-state temperature and grows logarithmically with system 
size. We investigate the effect of the other relevant parameters and estimate the force 
using the Fourier transform of the fluctuating hydrodynamic fields. 
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1. Introduction 

The Casimir force predicted in the seminal work of Casimir [1] is an attractive force 
F = —7T 2 hcA/ (240-D 4 ) between two perfect conducting neutral plates with area A facing 
each other at a distance D. This attraction originates from the modification of the long- 
range fluctuations of the quantum electromagnetic field, due to the boundary conditions 
imposed by the conducting plates [2j [3J [U [5]. The Casimir energy - the difference 
between the energies of the quantum electromagnetic field for the plates at distance D 
and at D—>oo, respectively - is proportional to h pQ HJ E] ; so is the Casimir force, that 
is the derivative of this energy with respect to the distance between the plates. 

Although the original Casimir interaction has a quantum nature, such an effect 
occurs in many classical systems where fluctuations are of thermal origin [2]. In fact, 
fluctuation-induced forces appear in systems with fluctuating long-range correlations 
that are geometrically confined by inserting external objects in the correlated medium. 
Examples can be found in nematic liquid crystals [6], critical mixtures (HE], superfluid 
films [9], and granular media [TOl [TT| fl"2] . In a thermally noisy correlated medium, 
where long range spatial correlations exist due to thermal fluctuations, the Casimir 
energy (and force) is expected to be proportional to UbT [2]. The results of a recent 
experimental study on the critical mixtures reveal that the Casimir energy is indeed 
very sensitive to the temperature of the system |13j . Depending on the characteristics 
of the system, the fluctuation-induced force has been found to be even repulsive, e.g. 
in dielectric materials with nontrivial magnetic susceptibility [14], or the interaction 
between a perfectly conducting and an infinitely permeable plate |15j . 

The fluctuation-induced force between two intruder objects immersed in a granular 
gas is studied in [10] , where it is found that the confinement of the fluctuation spectrum 
of the hydrodynamic pressure field induced by the intruders, leads to different local 
pressures in the gap between the intruders and the outside region. This effect causes 
an effective Casimir-like force between the intruders. The results of some experiments 
on granular fluids have confirmed the idea that the presence of large intruder particles 
modifies the thermodynamic properties such as pressure, velocity and density in the 
regions between the larger particles [T6l [T7J [18] . Appearance of long-range interactions 
(despite the short-range nature of the interactions on the grain scale) may shed new 
light on the mechanism of some collective behaviors, e.g. segregation [TTlfTB"] . in granular 
media. 

In this paper we present the results of extensive numerical simulations to study the 
fluctuation-induced force between two large objects immersed in a noisy granular gas. 
Although the interactions between particles in granular systems are dissipative [19j, we 
maintain the dynamics by means of an external driving force. Our main aims in the 
present work are to investigate the probability distribution of the fluctuating Casimir- 
like force and the relationship between the average force and important parameters of 
the system such as the steady-state temperature, the mass density of gas particles and 
the distance between the intruders. The question of whether the average force depends 
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Figure 1. Schematic picture showing the simulation cell. The pressure difference 
between the hatched and cross-hatched regions leads to an effective repulsive force 
between the two intruders. We choose L /r = 200, D /r — 30 and R /r = 10 in the 
reference system. 



on the system size is also addressed. 



2. Simulation Method 



The system we consider is a two-dimensional granular gas of identical rigid disks of 
radius r and mass m interacting via inelastic collisions. In order to exclude the undesired 
effects of side walls, periodic boundary conditions are applied in both directions of the 
square-shaped system of length L. Two immobile rigid intruder particles of radius R and 
infinite mass and moment of inertia are immersed in the granular gas bed, separated 
by a distance D. We have one reference system whose parameters are denoted by 
zero subscripts (see Fig. CD). Throughout the paper we either use this reference system 
or vary only one parameter to check its effect while other parameters are kept fixed 
at their reference values. We perform molecular dynamics simulations in which the 
number density n of the granular gas is 0.075r -2 and the effective normal coefficient of 
restitution a of the particle-particle and particle-intruder collisions is set to 0.8. 

The system is coupled to an external heat bath that homogeneously transfers energy 
into the system. The acceleration of each particle a, is perturbed instantaneously 
according to a[ = <Zj + where prime refers to the acceleration after perturbation, 
and £j can be considered as Gaussian white noise with zero mean and correlation 
(£ia(t)£jp(t')) = C,gSij5 a /35{t — t'), where a and f3 denote Cartesian components of vectors 
[2TH l2T] . Practically, the energy is transferred into the system in the following way 
[10], [20] : we update the momentum of each grain at each time step At. The components 
of the momentum are updated by adding random values that are chosen from a Gaussian 
distribution with zero mean and variance of. 
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3. Results 

The simulation results [solid line in Fig. [2(a)] reveal that the temperature T is time 
dependent and finally saturates. Note that the temperature T is a uniform field in 
our simulations in contrast to the case where energy flows into the system from the 
boundaries which leads to the spatial gradient of temperature [22]. We follow a mean- 
field approach to describe the time evolution of T. On the one hand the system gains 
energy due to coupling with the heat bath. The rate of the energy gain of a single 
particle averaged over the uncorrelated noise source g(£) is 

d t E = lim i- / [Ei(t + At) - ESMOdt = drn£. (1) 

At — >0 ZAt J £ 

On the other hand the system losses energy due to inelastic collisions. The rate of 
the energy loss of a sing le particle is described by d t E = -(1 - a 2 )ujT/d = (3T 3 ' 2 
[23] 124"] . where d is the dimension of the system, u((x yT) is the temperature- 
dependent collision frequency given by the Enskog theory [25J, and j3 is a coefficient 
that contains all of the relevant parameters except the temperature. Therefore the time 
dependence of temperature (energy) according to mean-field approximation is given by 
dT/dt = -f3T 3 / 2 + dm£ 2 . By integrating the equation dT/(-(3T^ 2 + dmQ = dt [26], 
one arrives at the following expression for the evolution of T [dash-dot line in Fig. [2](a)]: 

f(T)-f(T l ) = -dme o t/2T MF: (2) 

where Tj is the initial temperature and 



/(*) = l - In 1 - 2 ^- +3:/T - + ±= arctan ( 3 ) 

6 i + ^JtZ + x /t mf Vz 

Starting from any initial configuration, the driven granular gas finally converges 
to a nonequilibrium steady state, where energy dissipation due to inelastic collisions is 
balanced with the energy input. Although the mean- field prediction for the saturated 
temperature is T MF = (dm£ 2 /P) 2 ^ 3 , the temperature of the nonequilibrium steady 
state T NESS [dashed line in Fig. Eta)] is expected to be larger than T MF since it is 
logarithmically divergent in the system size (in 2D) due to the existence of spatial 
hydrodynamic fluctuations [20J. 

After the system achieves the stationary state, we study the effective interaction 
between the two fixed intruders. We measure the total momenta Pi and P r transferred 
from the granular gas to the left and right intruders respectively, during a time interval r. 
t corresponds to 10 4 time steps (~ 100 collisions per grain). The components of Pi and 
P r parallel to the line connecting the two centers of intruders provide the fluctuation- 
induced force as F=(P xr —P x i)/2r. Since this force is itself a fluctuating quantity, as 
expected for such correlated media [27], we let the simulation run until we measure the 
force for more than 10 4 consecutive time intervals r. Figure EJ^b) displays the fluctuating 
nature of F when measured in the steady state of the reference system. The probability 
distribution of F is plotted in Fig. Etc), where it turns out that the distribution can 
be well fitted to a Gaussian with the standard deviation 0.612T MF /r . The ensemble 
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Figure 2. (a) Temperature saturation starting from a random initial configuration 
of the system (solid line). The temperature T is scaled by the mean- field steady- 
state temperature T MF . The dash-dot line shows the time evolution of T according 
to Eq. Q. The nonequilibrium stationary temperature T NESS is shown with dashed 
line, (b) The time evolution of fluctuation-induced force F scaled by T MF /r a . (c) The 
probability distribution of F for the reference system where the simulation run over 10 8 
time steps. The dashed line shows the best fit with a Gaussian distribution centered 
on 0.015T mfV 



average of F in the reference system equals to (F) = 0.015T MF /r Q . The fluctuations are 
about two orders of magnitude larger than (F) , thus very long simulations are required 
to measure the force with small numerical errors. Such a repulsive force between two 
intruder objects sitting at a distance larger than their radii was reported in the previous 
studies of driven granular gases [HH [12] and even dense shaken granular packings [28] . 

The time step length - that reflects the time scale for interaction with the heat 
bath - is chosen large enough that the dissipation is kept alive, but it is much smaller 
than the mean free time between the collisions. Since the energy considerations yield 
£g = cr 2 /(m 2 At) [20], the steady state temperature (according to the Enskog theory) 
then becomes 

To investigate how the Casimir-like force is affected by the simulation parameters, we 
vary the parameter values one by one while the others are kept fixed at their reference 
values and measure the force. For this purpose, simulations are performed anew for 
each data point in Figs. E(a-d) and the force is measured after the system reaches the 
steady state. We note that the values of {F)/{F) around 5xl0~ 2 reflect the accuracy 
level of our calculations. Therefore, the best fits with power-law functions (not shown in 
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Figure 3. The fluctuation-induced force (F) scaled by (F) Q in terms of the grain 
radius r (a), the standard deviation of the momentum distribution a (b), the grain 
mass m (c), and the time step At (d). The dashed lines are power-law functions with 
exponents correspond to those of Eq. 0. (c) (F) versus the temperature T. Here, the 
dashed line displays a linear growth of (F) with T. 
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Fig. [3]) are obtained when the inaccurate data in Figs. [3](b) and (d) are neglected. Our 
results reveal that the exponents of the best fits approximately equal to the exponents 
of Eq. |4]with less than %6 errors. The stationary temperature is also measured for each 
data point in Figs. [3]^a-d). The average force is plotted as a function of the temperature 
in Fig. [3](e) that interestingly verifies the linear dependence of the fluctuation-induced 
force on the temperature, as theoretically was predicted in the literature for the thermal 
fluctuating correlated media [21 II]- 

Next we investigate how the effective interaction between the intruders is affected 
by the distance between them. Figure H](a) displays that the average repulsive force (F) 
decays with increasing the distance D. Supposing a power-law dependence of (F) on 
D, one gets the exponent —0.8 ± 0.3. 

The repulsive force originates from the pressure difference between the hatched and 
cross-hatched regions in Fig. [TJ The mechanism is briefly explained below (for details see 
[10J). In a system of hard disks, the pressure field p can be written as p(n, T) = TH{n) 
[29] , where n and T are fluctuating hydrodynamic fields. Expanding the pressure up to 
the second order around the stationary values (n ,T ) and taking its statistical average 
over the random noise source, we obtain (p) = p + Hi(Sn5T) +T MF if 2 (<5n 2 ), where 
p is the stationary pressure and Hi and H 2 are the first and second derivatives of 
H(n) with respect to n around no, respectively. By employing the Fourier transform 
of the fluctuating fields 8A(r) = J2 k e~ tkr 5A k /V and the structure factors Sab(^) — 
(5A k 5B^/V [20\ we rewrite p $0\: 

(p) = Pv + V-^WiSmity + T^HzSUk)}- (5) 
fc 

In the calculation of the pressure the main contribution comes from the region of small 
k. Substituting the structure factors of this region, Sab(}z) = S\ B /k 2 (where S AB is 
a function of the number density, restitution coefficient and noise intensity), we obtain 
the steady-state pressure: < p >= p +CV~ l J2k(^-/k 2 )i where C = HiS^ lT J r T MF H2S^ n 
is a function of the number density n and is negative for our simulations. In the gap 
between the intruders, the number of valid k modes (and therefore Ylk V^ 2 ) i s smaller 
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Figure 4. (a) The average force (F) scaled by (F) in terms of the distance between 
the intruders D scaled byr . The dashed (dash-dot) line corresponds to the force 
calculation with (without) taking into account the actual boundary conditions, (b) 
(F)/(F) as a function of the system size L. The dashed line corresponds to h(L) scaled 
by (F) . The inset shows how the deviation of the data from the logarithmic growth 
h(L) varies with the system size. 



than the outside region. Consequently, the pressure in the cross-hatched region is higher 
than the hatched region. This pressure difference Ap causes an effective repulsive force 
between the intruders. Interestingly, C (and therefore Ap and the fluctuation-induced 
force) is proportional to T MF . 

Ap is estimated in Ref. [ID] , supposing that k vectors are confined in rectangular 
boxes of size (D — 2R) x 2R and (L — D — 2R) x 2R in the gap between the intruders 
and in the outside region, respectively. The dash-dot line in Fig. 11(a) displays the 
Casimir force for different distances between the intruders, calculated according to the 
above mentioned estimation. The force is overestimated because (i) the geometrical 
simplification (using rectangular instead of circular boundary conditions) causes an error 
in the computed pressure, and (ii) the fluctuations are indeed correlated in the hatched 
and cross-hatched regions. We have improved this estimation by taking into account the 
circular shape of the intruders. At each point in the hatched or cross-hatched regions, 
the components of the allowed k vectors in the x and y (perpendicular to x) directions 
are {2,Tm x / djjj), 2irn y /L), where djjj)is the length of the line segment that horizontally 
connects the surfaces of the two intruders. Using the continuous form of 1/k 2 , the 
pressure at this point reads 

(2ir) 2 r 27T/r o r 2jT/r o i 

(p) = P +CV~^—L dkj dk y -—, (6) 

na x J2-rr/d x J2tt/L k x "+" K y 
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Where an upper cutoff 27r/r is used beyond which hydrodynamics is not valid. This 
cutoff is also needed to exclude the singular vector k = (0,0) from the integral. 
Calculating the pressure in the hatched and cross-hatched regions at a given y [according 
to Eq. (jSJ)], we obtain Ap(y). Next, by integrating Ap{y) over the whole range of y 
(-R, R) to obtain the average force, and repeating the process for different distances 
between the intruders we get the dashed line in Fig. IU(a), revealing a significant 
improvement due to taking into account the actual boundary conditions. We note 
that the correlation between the fluctuations in the two regions seems to be important 
since our results still overestimate the force. One expects that the interaction between 
intruders would be affected if a new intruder is inserted into the system. Therefore, our 
results provides an important evidence for the nonlinearity of the Casimir-like forces in 
thermally noisy systems and offers new insight into the problem which should trigger 
new experiments and theoretical explanations. 

Finally we study the system-size dependence of the effective force. It has been found 
in Ref. [20] that Tness—Tmf in a two-dimensional driven granular gas is logarithmically 
divergent in the system size due to hydrodynamic fluctuations. Such a size dependence 
is also expected for (F) since we have verified a linear dependence of (F) on T^ess- 
The simulation results [Fig. H](b)] reveal that (F) increases slightly as the system size 
L increases and the data can be well fitted to a logarithmic curve h{L) = a*lnL + b* 
with less than %4 error. Within the investigated range of L, the deviation of the data 
from the logarithmic relation does not show a systematic dependence on the system size 
[inset of Fig. H(b)]. The increase of (F) with L can be also explained qualitatively in the 
following way: When the system size is increased at a fixed D, the number of possible k 
modes (and 1/k 2 ) in the hatched region of Fig. [T] increases; consequently, the mean 
pressure in the hatched region is reduced while the pressure in the cross-hatched region 
is not varied. This leads to the increase of pressure difference between the hatched and 
cross-hatched regions which yields the increase of (F). 

4. Conclusions 

In summary, we have carried out numerical simulations of the fluctuation- induced force 
between two intruders in a driven granular gas bed. Earlier works have demonstrated 
that such a force is expected since the thermodynamic properties are affected due to the 
presence of the intruders. Here, our main focus has been to study the fluctuating nature 
of this interaction and the temperature dependence of the ensemble average of the force 
(F). We have verified that (F) increases linearly with temperature, and slightly with 
system size that can be described by a logarithmic growth. We have also improved the 
estimation given in [10J to explain the force, by taking into account the actual boundary 
conditions of the problem. 
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